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SUMMARY 


Modifications were made to the three-dimensional Douglas 
Neumann surface singularity method developed by Hess in order to 
improve the accuracy and efficiency of calculating the low speed 
potential flow on arbitrary lifting configurations. The original 
source paneling is replaced by a combination of source and doublet 
panels based on the classical Green's identity. The amount of 
calculation per panel is not increased, as the source distribu- 
tion is given directly by the body geometry. Solution to a set 
of equations for approximately one doublet strength parameter per 
panel then gives the doublet distribution to satisfy boundary 
conditions of flow tangency on the body. A doublet sheet repre- 
sents the wake of vorticity from the trailing edge of lifting 
surfaces . 

A numerical study of the characteristics of Green's identity 
was made for two-dimensional flow. The source panel method used 
in the original Douglas Neumann program is less accurate on thin, 
highly loaded lifting surfaces as the result of local source sin- 
gularities of large magnitude. Results show significantly better 
accuracy in such cases for the milder source-doublet combination 
of Green's identity. A wake-tangency Kutta condition improves 
the accuracy of the calculated lift. This permits a reduced 
number of panels (20 to 40 on an airfoil section) for a given 
accuracy, and a siibstantial saving of computing cost on config- 
urations such as high lift devices or supercritical airfoils. 

Other advantages of the Green's identity method are improved 
accuracy for the flow in concave corners and more reliable and 
efficient inverse solutions for design of multi-element sections 
with given pressure distributions. 

A comparison of different numerical methods for the Green's 
identity formulation resulted in selection of the approach used 
in the 3-D program modification. Flat, trapezoidal panels with 
piecewise constant source density cover the body surface and 
panels with piecewise quadratic doublet density cover the body 
and the wake. Boundary conditions of zero perturbation potential 
are satisfied at an internal control point on each body panel and 



at points on the edge of paneled sections. Special attention was 
given to reducing the discontinuity of doublet strength between 
adjacent panels. The internal doublet sheets used in the original 
program to provide lift carry-over through "non-lifting" bodies 
such as fuselages are not required in the modified program. The 
doublet distribution covers the outer surface of all body com- 
ponents, and the only distinction made between lifting and non- 
lifting components is in the joining of a wake to sharp trailing 
edges. 

The program is written for the CDC CYBER 175 at Langley 
Research Center. Calculated results are presented showing the 
accuracy and stability of the modified program for isolated 
bodies, wings, wing-body combinations, and internal flow. 
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INTRODUCTION 

The surface panel method philosophy for solvincr arbitrary 
incompressible potential flow problems involves the mating of 
classical potential theory with contemporary numerical techniques. 
Classical theory is used to reduce an arbitrary flow problem to 
a surface integral equation relating boundary conditions to an 
unknown singularity distribution (Reference 1) . The contemporary 
numerical techniques are then used to calculate an approximate 
solution to the integral equation. This involves representing 
flow boundaries by surface panels on which potential flow 
singularities are distributed. 

All properly formulated surf ace . panel methods are exact in 
the sense that the difference between the approximate numerical 
solution and the exact solution to the integral equation can be 
made arbitrarily small at the expense of increasing the number of 
computations. This does not imply that all panel methods are 
equally successful. Indeed, vast differences exist with respect 
to prediction accuracy versus computational effort, reliability, 
simplicity, and applicability to an inverse solution mode for 
design problems. 

The Douglas Neumann surface source singularity computer 
program (reference 2) has been used extensively at industry, 
university, and government facilities around the world for pre- 
dicting siobsonic inviscid flow. For most geometries, the program 
predictions are accurate, efficient, and reliable. However, for 
wings with thin trailing edge regions and high loading, the 
program predictions can be inaccurate or numerically unstable. 

The instability results from local source singularities of 
extremely large magnitude. 

This report describes the modification to the Douglas 
Neumann program to eliminate the difficulties associated with 
thin, highly loaded geometries. The basic concept involves 
replacing the source singularity formulation by the mild combined 
source-doublet distribution of the classical Green's identity. 
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The body of this report covers the theoretical fundamentals, 
the results of a study into numerical solution approaches for 
two-dimensional flow, the three-dimensional solution formulation, 
and representative numerical solutions for three-dimensional con- 
figurations. The section entitled Potential Flow Theory concen- 
trates on the concepts of importance in establishing a reliable 
solution formulation. Included in the section Research on Green's 
Identity Formulation are the effects in two-dimensional flow of 
(1) the type of surface singularity distribution employed, (2) 
higher order panel curvature corrections , and (3) the method of 
applying boundary conditions. Of special concern were predic- 
tion accuracy, numerical stability, efficiency, and applica- 
bility to inverse design problems. Using the two-dimensional 
study as a guide, the modification formulation for the arbitrary 
lifting body Douglas Neumann program was developed and is pre- 
sented in the section Numerical Solution Formulation . The 
example solutions in the section entitled Calculated Results 
include isolated bodies, wings, wing-body combinations, and 
internal flow. 
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POTENTIAL FLOW THEORY 

The solution to incompressible potential flow problems is 
based on the concept that arbitrary potential flov7 fields can be 
considered to be induced by suitable surface density distributions 
of singularities such as sources and doublets on the boundary 
contours. The utility of this concept is that simple flow fields 
induced by isolated singularities can be superimposed mathe- 
matically to generate fields having arbitrarily complex boundary 
geometries or boundary conditions. 

This section provides background on the classical funda- 
mentals of the theory, particularly the aspects pertinent to the 
formulation of a reliable and efficient numerical solution method. 
It will be explained that there is no limit to the number of dif- 
ferent surface source-doublet distributions that will induce any 
given flow field and that there is no theoretical basis to favor 
one distribution over another. However, such is not the case 
when the theory is coupled with contemporary numerical techniques 
to establish a solution method. The type of singularity distri- 
bution selected can play an important role in solution efficiency 
and numerical stability. An understanding of the theoretical 
fundamentals is necessary before it is possible to determine how 
best to model the theory in a numerical formulation. In that 
regard, the rationale for selecting the combined source-doublet 
distribution of Green's identity will be presented below. Sever- 
al two-dimensional examples that quantitatively demonstrate the 
advantages will be presented in the following section. 

Laplace Equation 

For the case of an incompressible, inviscid, and irrota- 
tional fluid, the Navier-Stokes equations can be reduced to the 
classical Laplace equation (reference 1) . Consider the fluid 
in the region R of figure 1. The illustrated fluid boundaries 
can be solid such as an airplane wing, permeable such as the sur- 
face of a jet stream, or imaginary. Because the fluid is incom- 
pressible, continuity dictates that the divergence of velocity 
must be zero. 
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The raomentuin equation is 


^ -)■ 
DV _ Vp 
Dt p 


(2 


where ^ is a material derivative, i.e., it refers to the time 
rate of change of a fixed set of fluid particles. By taking ti 
curl of equation (2), it is possible to establish Kelvin's law 
that an initially irrotational , inviscid fluid will so remain. 

“V 

The mathematical consequence of zero vorticity (V x V = 0) 
is that there exists a function $(x,y,z) which has a gradient 
equal to the velocity at each field point in the region R. 




V 


(3) 


= V$ 

Here it is assumed that the flow is initially irrotational , such 
as a free stream approaching a fixed airplane. The substitution 
of equation (3) into equation (1) generates the governing Laplace 
differential equation for incompressible potential, flow. 

V* (V$) = = 0 (4) 

Equation (4) applies whether or not the flow is time dependent. 

In the special case of steady state, Bernoulli's equation is 
readily established by integrating equation (2) along a stream- 
line path s between any two points A and B. 

B B 

f DV f ^ 

P ) Dt • " J • ds 

A A 

For steady flovi7, the differential ds can be replaced, by Vdt. 

B -> 

P ( Dt ' " ^Pb ~ Pa^ 

/V 2 V 2 \ 

P y 2 “ ~2 ) ~ ~ Pa^ 


p 


A 



P 





(5) 


If the region R is unbounded and point A is allowed to move in- 
definitely upstream, the conventional Bernoulli equation relat- 
ing pressure and velocity is established. 
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In the determination of the velocity field for irrotational 
flov7, the application of equation (4) is equivalent to the simul- 
taneous application of equations (1) and (2) . The obvious ad- 
vantage to equation (4) is that there exists only one unknown 
scalar function. For the usual case of steady flow, the solu- 
tion to the velocity field is essentially the complete solution, 
since equations (5) or (6) provide the pressure directly from the 
local velocity. Forces and moments can then be generated by 
pressure integration. The essence of the solution method in- 
volves determining the potential function $. 

Singular Solutions 

A simple function which satisfies the Laplace equation is the 
potential due to a point source of strength m positioned at 

<=‘o'yo’^o>- 


$ 


s 


-m 

4Trr 


(7) 


where 


1 

r s [(x-x^)^ + (y-y^)^ + 

2 

With the exception of the point (x^,y^,z^) , the Laplacian (v ) 
of equation (7) is zero. Since equation (4) need apply only to 
points in the flowfield, equation (7) satisfies the governing 
flow equation as long as (x^,y^,z^) is on or within the non- 
fluid side of a flow boundary. Equation (7) is designated a 
particular solution. The term "source singularity" refers to 
the fact that the field induced by a source satisfies the Laplace 
equation everywhere except at the singular point (x^,y^,z^). 

Because the Laplace equation is linear, the sum of any 
number of particular solutions is also a solution. Equivalently, 
any distribution of sources on or v/ithin flow boundaries is a 
solution to Laplace's equation. 
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other particular solutions can be generated by differen- 
tiating equation (7) with respect to . For example, 

the doublet function defined below, satisfies the Laplace 

equation except at the singular point 


t = s . (!^ 1!£'1 ( 8 ) 

\»Xo' 3y^’ 

The direction of unit vector n is arbitrary. Physically, equation 
(8) describes the limiting situation of an equal strength source 
and sink on an axis parallel to n (figure 2). The separation dis- 
tance is equal to As and the source strength is 




Figure 2. Point Doublet at (xq, Vq, Zq) 
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Any distribution of sources and doiablets on or within flow 
boundaries automatically satisfies the governing Laplace equa- 
tion. It is possible to demonstrate that any boiinded potential 
flow field can be considered to be induced by a continuous dis- 
tribution of sources and doublets (reference 1) . Indeed, there 
is no limit to the number of different distributions that will 
induce any one given flow field. Once such a distribution is 
found, the solution is essentially complete because the flow 
potential and velocity at any point can be calculated directly 
as the integral sum of the individual differential contributions 
induced by the source and doublet distributions. 

This study is restricted to the investigation of source and 
doublet singularity distributions on flow boundaries, not within. 
An example of the latter is the representation of the flow in- 
duced by a solid body of revolution immersed in a uniform stream. 
For some body shapes, it is theoretically possible to represent 
the effect of the body by an axial distribution of source sin- 
gularities (figure 3). However, in the formulation of a numer- 
ical solution method, a distribution of singularities on the body 
surface provides a far more attractive model. It is desirable 
to have the singularities as close as possible to the flow which 
they induce. Otherwise, the distinction between the influences 
of neighboring regions of the singularity distribution is clouded, 
and numerical instabilities can result. This effect becomes more 
pronounced as body thickness increases. 



Figure 3. Body of Revolution Represented by an 
Equivalent Axial Singularity Distribution 
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Possibly a more compelling reason for avoiding internal 
singularities regards the existence of solutions. For bodies 
with flat disk noses or tails such as in figure 4, there is theo- 
retically no axial distribution of sources that can induce the 
flow. On the other hand, there always exists a solution 
surface distribution of singularities on the flow boundaries. 



Figure 4. Body for Which IMo Equivalent Axial 
Singularity Distribution Exists 

In accordance with equations (7) and (8) , it is possible to 
express the induced potential $ at any field point P as the in- 
tegral sum of continuous source and doublet distributions on 
boundary surfaces (see figure 5) . 

» = J {- * il- 

In equation (9), the following definitions are used: 

(1) a and \x are respectively the source and doublet 
strength per unit area on the boundary at arbitrary 
point Q, 

(2) dS is a differential surface area element at point Q, 

(3) r is the distance between boundary point Q and field 
point P, and 

(4) n is distance measured along an axis normal to the 
boundary surface at Q, positive into the flow. 

Notice that the doublet axis direction in equation (9) is normal 
to the flow boundary. 
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Figure 5. Points in a Fluid and on a Boundary 

There is an important equivalence between surface doublet 
distributions and vorticity. As is proved in reference 2, a 
surface doublet distribution of density y can be replaced 
by an equivalent surface vortex distribution where the vortex 
density vector y ^.t each surface point satisfies the following 
equation : 

Y=nxVy (10) 


where n is to be interpreted as the unit normal vector. Equa- 
tion (10) follows from the fact that a vortex loop of strength 
r induces the same velocity at every field point as a uniform 
doublet sheet of density r, provided that the edges of the 
sheet lie on the loop (figure 6). The theoretical equival- 
ency between vorticity and doublet distributions does not imply 
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equivalent simplicity in a niimerical formulation. In this study, 
attention will be focused on doublets rather than vorticity, 
because distributions of scalars such as y are generally much 
simpler to handle than vector distributions such as vorticity. 




Constant Density Doublet Surface 
(y ->• Strength/Area) 


Figure 6. Equivalency Between Vortex Loops and Doublet Surfaces 

For any bounded potential flow field, there exists a surface 
source distribution a and doublet distribution y that will in- 
duce the field. Integral equation (9) describes the resultant 
potential at any field point P, and in accordance with equa- 
tion (3) , the gradient of equation (9) will give the flow velo- 
city at point P. As discussed above, equation (9) satisfies 
the Laplace equation regardless of distributions a and y. It 
remains to determine the appropriate distributions, and this is 
accomplished by satisfaction of boundary conditions. 

Boundary Conditions 

Within certain constraints on geometric slope continuity 
(reference 1) , a bounded, simply connected velocity field is 
uniquely determined by the distribution on flow boundaries of 
either the normal component of total flow velocity V $ • n or 
the total potential $ . These boundary value problems are res- 
pectively designated Neumann and Dirichlet problems in the clas- 
sical literature. In this study, primary consideration will be 
given to Neumann problems, because most practical cases involve 
prescribed normal velocity boundary conditions. In particular. 
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the flow tangency associated with solid bodies is expressed 
mathematically as V$-n = 0. All formulations in this study 
are sufficiently general to allow direct extension to arbitrary 
Ne\amann, Dirichlet, or mixed Neumann-Dirichlet boundary condi- 
tions . 

The relationship between prescribed Neumann boundary condi- 
tions and the unknown source-doublet distribution is generated 
by taking the gradient of equation (9) and then allowing field 
point P to approach the boundary surface at arbitrary boundary 
point B. The normal velocity component at B is designated 






> ( 11 ) 


Vjq = lim V $ . ng 

P-^B ^ 

->■ 

where is the unit normal vector at boundary point B. 

The limiting procedure avoids contacting the singularities; this 
is necessary since the Laplace equation is satisfied arbitrarily 
close to a source or doublet location, but not at the location 
precisely. 

There exists a solution distribution o and p for any arbi- 
trary normal velocity distribution provided that § ~ 

The integral constraint expresses the fact that the net fluid 
mass surrounded by the boundaries cannot change with time for 
incompressible flow. At each boundary point B there are two 
unknowns, a and p, while there is one boundary condition, the 
prescribed value As one might expect, for any set of 
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boundary conditions, there is no limit to the number of dif- 
ferent solution combinations of source-doublet distributions. 

In other vvords, there is no limit to the number of different 
surface singularity distributions that can theoretically induce 
the exact same field. The discussion to follow is intended to 
establish the particular source-doublet combination most suitable 
for application to a numerical solution formulation. 

One possible candidate is a source -only distribution, with 
y set to zero a priori. The disadvantage is that source-only 
solution methods can generate strong numerical instabilities 
and unreliable predictions. The explanation is that for thin 
bodies at incidence, the solution source distribution tends to 
increase in magnitude linearly with respect to the inverse of 
thickness, while the net flow velocity changes only slightly. 
Prediction errors tend to be proportional to the singularity 
strengths. For a thin body, the average magnitude of source 
density can be several times that of the induced velocity, with 
a subsequent magnification of prediction error. An analogous 
situation could occur in m.easuring the distance AB in figure 7. 

The distance could be determined indirectly by subtracting the 
measured distance BC from measured AC. In the absence of mea- 
suring errors, AB is exactly AC - BC. However, if the relative 
error in taking a measurement is e, the indirect measurement would 
generate the following error AAB : 


A 


Direct: 


Measure AB 


Indirect: 


TMeasure AC and BC 




AB = AC - BC 


} 




Small Error 
Large Error 


Figure 7. Direct and Indirect Measurements 
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AAB ~ „ . AC 

= ze 

AB AB 

If AC is sufficiently large, the error in measuring ^ would be 
greater than 100%. Obviously, one would be far better off mea- 
suring the distance AB directly. In the same spirit, one would 
prefer to use singularity distributions of approximately the same 
magnitude as the induced flow velocities. 

This can be accomplished by applying Green's third identity 
(ref. 3, p.219). If the source distribution is established a 
priori by setting a equal to at each boundary point, and then 
if equation (11) is satisfied, it can be proved that m = at 
each boundary point. By making these substitutions for a and y 
in equation (9), Green's identity is established. 

* ^ 17 3^(|)}‘3S 

The utility of Green's identity is that the singularities are 
no stronger than the flow field which they induce; indeed, there 
is a local equality regardless of geometry or boundary conditions. 
This contrasts sharply with source only solutions for thin geome- 
tries subjected to high lift. For the flow around conventional 
slotted flap high lift devices, the source densities can be orders 
of magnitude greater than the induced flow velocities. 

The demonstration that setting a = results in y = $ is 
straightforward if consideration is given to the imaginary flow 
field inside the real flow boundaries (figure 8) . Subscripts 
E and I are used to designate the real external and imaginary 
internal fields, respectively. The following discontinuities 
apply to any sheet of sources and doublets. 


‘^E " ^I = ^ 


(13) 
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E ^ Fluid Side of Boundary 

(External) 



Figure 8. Two Sides of a Flow Boundary Surface 


Both VNg and are considered positive toward the real flow 

(direction n) . Since a = VN^f it follows from equation (13) 
that VjsIj = 0. The imaginary flow must be stagnant because the 
normal velocity component at all boundary points is zero and 
because such Neumann boundary conditions generate a unique solu- 
tion. Therefore, is a constant, which can be selected as zero. 
Finally, equation (13) indicates that $„ = y , completing the 

III 

demonstration . 
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In most problems of interest the flow is unbounded. In that 
case, the surrounding boundary is allowed to increase indefinitely 
and its effect is replaced by adding a free stream of uniform 
velocity . However, Green's identity, as expressed by equation 
(12), still applies. Indeed, if solid bodies are immersed in the 
free stream, equation (12) indicates that a = 0 on the solid body 
surfaces and a doublet only representation results. As will be 
demonstrated later, convergence difficulties arise when doublet 
methods are applied to design problems. 

By applying Green's identity to the perturbation potential 
instead of to the total potential 3>, a particularly useful com- 
bined source-doublet representation is defined (Reference 1) . The 
representation will be shown to provide reliable solutions for both 
analysis and design problems. The perturbation potential tj) is 
defined as follows (see figure 9) : 

$ = V^* (x,y , z) + 4> (14) 



Figure 9. Body Immersed in an Unbounded Flowfield 
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To apply Green's identity to the perturbation potential, 
is replaced by V(J>-n (rather than V$-n) and y is replaced by <|) 
(rather than $) . 


■#> = iT {- 


q , y 3 
4irr 4it 3n 



a 



(15) 


where V(j>.n = ^co is a simple function of 

local geometry. Because <j> approaches zero at infinity, there is 
no need to include a surrounding boundary, which considerably 
simplifies the formulation of a numerical solution method. In 
equation (15) the singularity densities a and y are equal to 
the perturbation normal velocity and perturbation potential, 
respectively. If the perturbation flow vanishes, so do the 
singularities which induce the flow. Hereinafter, the term 
"Green's identity" refers to equation (15). 

The singularities of Green's identity are mild, and there 
are no more unknowns than for a source only formulation. For 
Neumann problems the source density distribution is known a priori 

a = V(j).n = Vn - V,„.n 

and the doublets are to be determined. For Dirichlet problems, 
the opposite occurs; the doublet distribution is known from the 
prescribed potential distribution on the boundaries, and the 
source distribution is to be determined. The approach of apply- 
ing Green's identity to a numerical solution method is believed 
to have been first employed by Morino, et al (reference 4) . 

The equivalency between a doublet distributon and vorticity 
as expressed by equation (10) results in a useful equality for 
Green's identity. Since y = <() on boundary surfaces and since 
^ = ^oo + V(j) , it follows that the tangential component of pertur- 
bation velocity on the surface will be equal in magnitude and 
normal in direction to the vorticity density vector y* Figure 
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10 illustrates the situation for the special case of two-dimen- 
sional flow. The total flow velocity vector V is seen to con- 
sist of components and where 


Vt = + y 

Vn = %„ + a 


(17) 



An illustration of the nature of three different singularity 
distributions which induce exactly the same flow field is present- 
ed in figure 11. The problem is to calculate the flow about a 
solid two-dimensional circular cylinder immersed in a uniform 
stream . The analytical singularity distributions are pre- 
sented for a source only solution, doublet only solution 
(Green's identity applied to total potential $) , and the combined 
source-doublet distribution corresponding to Green's identity 
(applied to perturbation potential 4)). For clarity, the equival- 
ent vortex density y is presented instead of doublet density. 
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It is noteworthy that the combined source-vortex solution is 
the average of the source only and vortex only solutions. 
Green's identity typically provides a source distribution more 
mild than the source only solution and a vortex distribution 
more mild than the vortex only solution. 


Representation 


Source distribution only 
— — ^ ”(B) Vortex distribution only 
' ■ — —“(6) Green's identity 


SOURCE 

DENSITY 

a 


VORTEX 

DENSITY 

7 
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Figure 11. Equivalent Singularity Representations 
for a Circular Cylinder 

Green's identity also offers the prospect for using internal 
indirect boundary conditions (reference 4) . In solving Neumann 
problems, the first step is to set a = Vjq - ^„.n, where 
* n is a function of geometry only. The second step for 
direct boundary condition satisfaction is to solve for y by 
satisfying equation (11) at each boundary point. However, it 
can be shown to be theoretically equivalent to establish zero 
perturbation potential on the imaginary or internal side of all 
flow boundaries. The proof follows immediately from equation 
(13) . 
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The advantages -of Green's identity in the application to a 
numerical solution method are summarized below. 

(1) Mild combined source-doiiblet distributions suppress 
numerical instabilities, which can otherwise be pre- 
valent for thin high lift geometries. 

(2) Singularity distributions vanish as the perturbation 
field vanishes, thereby eliminating possible residual 
error. 

(3) Direct relationships between velocity and singularity 
strengths on boundary surfaces exist which simplify 
calculations . 

(4) Possibility exists of applying indirect internal 
boundary conditions to simplify calculations. 

One additional advantage is in the application to design or 
inverse problems, for which the geometry most nearly correspond- 
ing to a prescribed velocity distribution is to be determined. 

Even though doublet only analysis problem solutions are rela- 
tively mild, they tend to be unsuitable for inverse problem 
application in leading edge regions. Supporting two-dimensional 
examples and the explanation will be furnished in the next section. 

Wakes 

It is possible to prove that the net force on all bodies 
immersed in a steady, incompressible, potential flow field is 
zero (reference 1) . In other words, an airplane would have no 
lift in a fully potential fluid. On the other hand, it is common- 
place to solve aerodynamic problems by representing lifting 
airplanes as solid bodies immersed in a potential fluid. The 
representation is made reasonable by the introduction of an 
additional flow boundary, the semi-infinite trailing wake. 

The need for a wake can be illustrated by considering 
flow circulation. Consider a section of an airplane wing 
(figure 12). The circulation r is defined as 

r E ^ V.dl (18) 

A" 
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r 



Figure 12. Circulation Around a Wing Section 

where the path of integration is a closed contour in the fluid 
surrounding the wing. If the fluid were to satisfy Laplace's 
equation everywhere, then equation (3) could be substituted into 
equation ( 18 ) . 

r = ^ V$.dl = 3)^+ - 

The circulation would be zero because <I> would have to be con- 
tinuous at A to s^atisfy Laplace's equation. Of course, a lifting 
wing possesses circulation, and therefore the fluid cannot 
realistically be described by Laplace's equation alone. The 
missing ingredient is the effect of viscosity at sharp edges. 

Potential fluids theoretically generate infinite velocity 
gradients at sharp convex corners such as at wing trailing edges. 
Regardless of free stream Reynolds number, the action of vis- 
cosity will be felt if velocity gradients are sufficiently 
s'trorig. Therefore, even what is conventionally designated an 
inviscid fluid will demonstrate the effects of friction at a 
sharp trailing edge. The fluid will search for a finite velocity 
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at the sharp edge, corresponding to what is referred to as 
the Kutta condition. The fluid in this localized trailing edge 
viscous region will be carried downstream to form a region v;hich 
cannot be described by Laplace's equation. This viscous region 
is, of course, the wake. 

As free stream Reynolds number increases, the thickness of 
the wake region diminishes and is often assumed to be zero. 
Nonetheless, the fluid properties are discontinuous across the 
wake and the wake itself cannot be neglected. Usually the wake 
is modelled as a stream surface having zero pressure loading 
imposed by the surrounding flow. These two free stream surface 
conditions are sufficient to determine the two sets of unknowns, 
wake geometry and wake singularity distribution. 

A frequently used assumption for wake geometry that is 
approximately correct for sufficiently high Reynolds numbers and 
small flow perturbations is that the wake is a flat surface 
emanating from the trailing edge and parallel to the free stream. 
The condition of zero loading is represented by the small distur- 
bance approximation of zero vortex strength along lines normal to 
the free stream (Figure 13). Equivalently, the streamwise gradi- 
ent of wake doublet density is assigned the value zero. 

The two-dimensional flow analogue to a wake is a streamline 
along which velocity is continuous but potential is discontin- 
uous. The strength of the discontinuity is equal to the circula- 
ation and is dictated by the Kutta condition of finite trailing 
edge velocity. In its simplest form, the Kutta condition is to 
be considered an approximation to the viscous part of the flow 
problem. 

In this study, the wakes are represented by a doublet sheet 
having zero thickness and zero streamwise doublet gradient. The 
wake geometry, though not necessarily flat, is specified. In the 
case of multi-component wings such as slats and slotted flaps, 
each component is assigned a separate wake. 
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RESEARCH ON GREEN'S IDENTITY FORMULATION 


A study was made in a two-dimensional surface paneling 
formulation to give a better understanding of the Green's 
identity method, as a basis for selecting the numerical methods 
to be used in applying the Green's identity formulation in the 
3-D program. A large number of formulations could be compared 
using available programs, and essentially "exact" calculations 
could be made by transformation methods (references 5, 6) to 
compare the accuracy of the different methods. 

The classical potential theory of the previous section 
reduces an arbitrary flow problem to a surface integral equation 
relating boundary conditions to an unknown singularity distri- 
bution. Numerical techniques are then used to calculate an 
approximate solution to the integral equation. Paneling methods 
do this by dividing the surface into a number of panels with a 
source or doublet singularity density distribution of unknown 
magnitude on each panel. The singularity densities are then 
calculated by solving a set of equations satisfying boundary 
conditions at control points on each panel plus additional Kutta 
conditions for the circillation to give smooth flow at Kutta 
points on the trailing edge of lifting surfaces. 

Many conceivable approaches are available for formulating 
surface paneling methods (references 2, 4, 7-18). At the level 
of potential theory, an unlimited number of combinations of 
source and doublet singularity distributions can induce a given 
external flow field. Our main interest in this study is in the 
particular source-doublet combination defined by Green's iden- 
tity, as well as the source only method used in the original 3-D 
Douglas Neumann program. Additional choices can be made of the 
numerical techniques to be used. The results of a study into 
the merits of the various approaches are presented in this 
section . 
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Comparison of Surface Paneling Methods 

The principal distinguishing features of surface paneling 
methods are (1) choice of representative singularity distribu- 
tion, (2) numerical scheme for distributing the singularities, 

(3) type of boundary conditions employed, and (4) conversion of 
solution singularity densities to velocity. A study was made 
to determine which selections in the above four categories pro- 
vide the most desirable characteristics in terms of prediction 
accuracy, computational efficiency, and numerical stability for 
both analysis and design. 

Special attention was given to problems for which the Green's 
identity formulation showed promise of signficant improvement 
compared to the original Douglas Neumann source method. The 
Douglas Neumann method uses flat, constant density source panels, 
plus a constant vorticity (or linear doublet) distribution 
around an airfoil or wing to give the circulation which satis- 
fies a Kutta condition of equal pressure on the upper and lower 
surface panels next to the trailing edge. Reference 8 discusses 
some problems with this method resulting from, large source den- 
sities on thin, highly loaded surfaces. The accuracy can also be 
poor when the panel number around a section is reduced to the 
range of 20 to 40, typical of 3-D applications, due to inaccurate 
circulation and panel curvature effects. As discussed later in 
this section, the Douglas Neumann- method erroneously predicts un- 
bounded velocities as a sharp, concave corner is approached. 

These characteristics present difficulties in applying the method 
to problems of current interest, such as supercritical airfoils 
and strongly cambered high lift devices such as slats and vanes. 
Accordingly, the study included analysis and test cases related 
to (1) accuracy vs panel number and panel distribution, (2) 
thin surfaces, with a small value of the ratio of distance 
between panels to panel length, (3) method of applying the Kutta 
condition, and (4) flow in concave corners. 
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Table I indicates the combinations of formulation approaches 
which were tested on two-dimensional sample problems during the 
study. The designations "low order" and "high order" refer to 
the numerical scheme for modeling the surface singularity dis- 
tributions- In every case, the low order scheme incorporates 
a piecewise constant strength source distribution, a continuous 
piecewise linear vortex strength distribution, and flat surface 
panels. The high order scheme includes corrections for local 
source gradient and panel curvature in accordance with reference 
8. Boundary conditions are normally applied directly at one 
control point per panel. For both the doublet and Green's 
identity approaches, theoretical relationships exist between 
local singularity strength and local velocity which can be used 
to provide either an equivalent internal potential boundary 
condition in lieu of prescribed normal velocities or a convenient 
means for calculation of velocity from singularity strength. 


TABLE I. 2-D SURFACE SINGULARITY METHODS TESTED 



Surface Paneling 

Boundary Conditions 

Velocity Calculation 


Singularity 

Low 

Order 

High 

Order 

Normal 

Velocity 

Internal 

Potential 

Influence 

Coefficient 

Summation 

Local 

Singularity 

Strength 

Inverse 

Solution 

Mode 

1. Source 

X 

X 

X 


X 



2. Doublet 

X 


X 


X 

X 

X 

3. Combined 
Source-Doublet 


a. Least Squares 
(Reference 17) 

X 


X 


X 


X 

b. Symmetric 
Singularities 
(Reference 18) 

X 


X 


X 


X 

(1) 

c. Green's 

Identity (2) 

X 


X 


X 

X 



X 

X 


X 

X 


(Ref 19-21) (3) 

X 



X 


X 

X 
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T;vpe or oounaar 

Existing computer programs were available to MCAIR for all 
the combinations listed in Table I, except for the application 
of internal potential boundary conditions to a Green's identity 
formulation (which was developed as part of this contract effort) , 
and the inverse (or design) capability for the Green's identity 
formulation with internal potential boundary conditions (which 
was recently developed in a separate MCAIR supported effort.) 

Documentation of the earlier programs is available in references 
17-20, as noted in Table I, and the two new Green's identity 
programs are discussed later in this section (also see Reference 21) . 

Detailed discussion and supporting examples comparing most 
of the formulation approaches indicated by Table I have been 
provided in references 17, 19, and 20 and will not be duplicated 
here. However, the more significant characteristics are repeated 
below. 

1. Source methods are numerically unstable for thin highly 
loaded geometries, which results from using unnecessarily 
strong singularity densities to induce a flow field. 

(See figure 14) . 

2. Doublet methods fail near leading edges in design 
problem iterative solutions. This is the consequence 
of the theoretical equivalence between doublet gra- 
dient and local velocity. In low speed regions, the 
combination of low flow velocity and low singularity 

strength numerically decouples the effect of local ' 

geometry from velocity, and the result is numerical 
instability. (See figure 15). 

3. The combined source-doublet combination obtained by 
the least squares approach of reference 17 is nximeri- 
cally stable for both analysis and design. However, 
the least squares procedure results in a significant 
increase in computational expense by increasing the 
n\amber of unknowns. 

4. The symmetrical singularity approach of reference 18 
is unstable for design. 
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Figure 14. Singularity Distributions for Supercritical Airfoil at 4° « 


30 






5 . 


The mild combined sour ce-doiob let distribution corres- 
ponding to Green's identity generates numerically 
stable results with none of the additional expense 
associated with least squares. 

6. The high order niimerical scheme applied to flow tan- 
gency boundary conditions provides increased accuracy 
compared to low order at only a slight increase in com- 
putational expense and with no increase in numerical 
instability. (See figure 16) . 

7. If flow tangency boundary conditions are applied to 
the high order Green's identity method, the velocity 
calculated by summing the individual contributions of 
each panel is significantly more accurate than applying 
the relationship between local velocity and local 
singularity density. 



Figure 16. Two-Element Airfoil Inviscid Solution 
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8 . 


The mild singularity distribution near the trailing 
edge from Green's identity permits improved accuracy of 
circulation, by a Kutta condition of wake tangency (zero 
normal velocity) at a point on the trailing edge bisector 
a few percent of the trailing edge panel length down- 
stream of the trailing edge. 

Internal Potential Boundary Condition Formulation 
Inasmuch as no computer program which could handle internal 
potential boundary conditions was available to MCAIR at the 
start of this contract study, the above enumerated characteris- 
tics of surface singularity methods do not include consideration 
of potential boundary conditions. Therefore, an appropriate 
computer program was formulated and tested during this reporting 
period. The remainder of this section is dedicated to a brief 
explanation of the program formulation, the presentation of 
illustrative examples, and the resulting conclusions. 

As implied by equation (17) and illustrated in figure 10, 
Green's identity implies that if the source density a everywhere 
on a geometry satisfies the relationship 


a 




N 


(19) 


then it follows that the vortex density y will satisfy the fol- 
lowing equation 


T = Vt - Vt (20) 

CO 

Here subscripts N and T refer to the surface normal and tangen- 
tial directions respectively, V is velocity, and refers to free 
stream conditions. It is simple to demonstrate that if equations 
(19) and (20) are satisfied, the flow field inside the body is 
uniform with velocity V = (reference 19) . Equivalently, the 
disturbance potential <}) inside the body will be a constant, and 
this fact can be utilized to generate equivalent internal poten- 
tial boundary conditions. The procedure, first applied by Morino, 
et al, (reference 4) is as follows: 
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1. Calculate the source density a on each panel from 

equation (19) . is known a priori from prescribed 

Ne\amann boundary conditions and is a function of 

local surface slope only. 

2. Calculate the local vortex density on each panel by 
solving a system of linear equations corresponding to a 
constant internal perturbation potential at each 
panel control point. 

3 . Calculate the velocity at any field point directly from 
the solution source-vortex distribution. 

For airfoils having sufficiently small camber and thickness, 
it can be demonstrated that the above solution procedure is equiv- 
alent to small disturbance,, linearized boundary condition flow 
theory. The demonstration is based on the fact that along any 
internal path connecting any two control points the average pertur- 
bation velocity component parallel to the path is zero. This is 
a direct consequence of having the same perturbation potential at 
all internal control points. Along a path connecting upper and 
lower surface control points at a given chordwise station, the 
average perturbation velocity component normal to the chord will 
be zero. The source contribution to this zero average is approx- 
imately equal to the free stream velocity component normal to the 
local camber, a fact which is easily verified by considering the 
difference between the upper and lower surface source densities 
(Equation 19) . The source induced normal velocity must be offset 
by the vortex contribution in order to attain the magnitude of zero. 
Therefore the average of the upper and lower surface vorticity is 
uniquely determined by the camber and incidence in the same manner 
as for linearized flow theory. Thickness effects can be understood 
by considering the average perturbation velocity along a chord- 
wise path between two internal control points on either the 
upper or lower surface. Because the average of the upper and 
lower surface sources is approximately equal to the thickness slope 
(Equation 19) , the source contribution to the internal chordwise 
velocity is the same as the thickness-induced tangential velocity 
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of linearized flow theory. The vortex- induced internal velocity 
along the chord line will offset the source contribution, with the 
result that the difference between upper and lower surface 
vorticity is twice the thickness-induced velocity of linearized 
flow theory. The consequence of the similarity between the appli- 
cation of internal potential boundary conditions and small 
disturbance linearized flow theory is that for sufficiently thin, 
planar geometries, both approaches will generate similar numerical 
behavior. In that regard, the uniform strength source-doublet 
quadrilateral panel method of Morino (Ref 4) actually reduces to 
a conventional vortex lattice method as thickness approaches zero. 



Figure 17. Two-Dimensional Panel Modeling 

(Flat Panels) 

In the present formulation, a low order modeling was used 
(piecewise uniform source, continuous piecewise linear vorticity, 
flat panels) . See figure 17. At the center or control point of 
each panel i, the internal perturbation potential is designated 

Uniform internal perturbation potential at each control point 
is obtained by satisfying the following system of equations: 

^<(>1 = = 0 (l£i£n-l) (21) 
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where n is the number of panels used to model the airfoil element. 

An additional equation is introduced to fix the circulation. At 
the discretion of the user, the circulation can either be prescribed 
explicitly or generated implicitly by a Kutta condition. The un- 
knowns are the vortex densities at the n panel endpoints. The. 
source densities are known a priori from equation (19) . 

After solving the system of linear equations for the unknown 
vortex densities, the surface velocity can be computed either by 
summing the effects of the individual panel singularity distribu- 
tions or by resorting to equations (19) and (20) , which provide 
a direct relationship between local velocity and local singularity 
density . 

In order to solve equation (21) ,■ it is first necessary to 
establish the linear relationship between potential and the source- 
vortex distributions on a panel. Consider a panel of length s in 
Figure 18. The following equation provides the potential at an arbi- 
trary point ( induced separately by a constant source density, 

constant vortex density, and linear vortex density on the panel. 



Figure 18. Coordinate System for a Panel of Length s 
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It is noteworthy that the vortex induced potential is dis- 
continuous at tIq = 0 for 5 ^ > 0. In establishing Acf)^ for equation 
(21) , it is necessary to eliminate the effect of the discontinuity. 
This can be accomplished by interpreting A<j)^ as the line integral 
of velocity along any internal path between control points i and 
i+1. For convenience, the selected path follows the surface of 
panels i and i+1. If this path crosses Rq = 0 for 5 ^ > 0, then 
a value equal in magnitude to the potential discontinuity is 
subtracted to eliminate the effect of the discontinuity. 
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Special treatment is applied if there are any slope discon- 
tinuities on the airfoil surface, such as a sharp trailing edge. 

At each panel endpoint that is designated a sharp corner, a dis- 
continuity is allowed in the vortex density. This corresponds 
to the fact that theoretically the perturbation velocity:' is dis- 
continuous at a corner. The magnitude of the vortex density hop 
at the panel endpoint is an additional unknown which is determined 
through the introduction of an internal control point a few percent 
of panel length from the corner. The imposition of uniform 
internal perturbation potential is applied to each sharp corner 
control point, as well as to the panel center points. 

Computing time is nearly independent of the type of boundary 
condition applied. For either the direct flow tangency or 
indirect internal potential methods, essentially the same terms 
are required for influence coefficient calculation and the same 
number of linear equations must be solved. 

Comparison of Green's Identity Formulations 

The present low order Green's identity panel method using 
potential boundary conditions has been compared with both the low 
and high order Green's Identity formulations of references 19 and 
20, in which the conventional prescribed normal velocity boundary 
condition is applied at each panel control point. It is noted 
that in every test case for which velocity or pressure is 
illustrated in the results to be presented below, equation (20) 
was applied for the potential boundary condition method and the 
individual effects of the panels were summed for the other two 
methods. This corresponds to what has proved to generate the 
best prediction accuracy in each case. In selected cases, re- 
sults are also presented for the low order Douglas Neumann 
source method which uses piecewise constant source density, 
uniform surface vorticity, and flat panels. This source method 
is analogous to the existing three-dimensional Douglas Neumann 
Program of reference 2. In every example for which there is 
lift, the same net integrated vorticity was used for all the 
methods in order to avoid the clouding effects of different 
Kutta conditions . 
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The test case solid-body geometries cover a wide range of 
shapes and thicknesses. Included are a circular cylinder, 
Karman-Trefftz airfoil, supercritical airfoil, nearly flat 
plate with noisy coordinate definition, and a thin symmetrical 
airfoil at zero incidence. The symbols used in the figures for 
all cases are defined as follows: 

O Low order Green's identity with potential boundary 
conditions 

A Low order Green's identity with normal velocity boun- 
dary conditions 

□ High order Green's identity with normal velocity 
boundary conditions 
<> Low order Douglas Neumann 

For the example involving a circular cylinder, the objective 
was to determine prediction accuracy versus panel density. Two 
panel models were established, one using twenty and the other 
using forty equally spaced panels. By comparing the calculated 
surface vortex density y against the exact analytical value 
corresponding to Green's identity, it was possible to establish 
the root-mean-square error in vorticity around the cylinder. 

The results are as anticipated for the case involving prescribed 
normal velocity boundary conditions, i.e., flow tangency (figure 
19). Consistent with the explanation of reference 20, the low 
order prediction error is proportional to the inverse of panel 
density and the high order error is proportional to the square 
of the inverse of panel density. What was unforeseen (as shown 
in figure 19) is that the prediction error for the low order 
solution with internal potential boundary conditions exhibits 
higher order characteristics; that is, the error varies propor- 
tionately to the square of the inverse of the number of panels. 

Similar trends are apparent for a thin, highly loaded 
Karman-Trefftz airfoil (figures 20 and 21). The virtually 
exact conformal mapping technique of Catherall-Sells (reference 
5) was used to provide a standard of comparison. Both in the 25 
panel representation of figure 20 and in the stagnation region 
of a 50 panel representation (figure 21), the calculated velocity 
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O Potential boundary conditions (low order) 
A Flow tangency (low order) 

□ Flow tangency (high order) 



Figure 19. Prediction Accuracy for a Circular Cylinder 


and pressure distributions for the low order potential boundary 
condition solution are very accurate and nearly identical to the 
high order flow tangency solution. On the other hand, the 
accuracy of the low order flow tangency solution is considerably 
worse . 
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Figure 20. Karman - Trefftz Airfoil 
0° a, 3.24 Cg, 25 Panels 
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For the supercritical airfoil of Figure 22, the low order 
potential boundary condition method again provides excellent results. 
The Douglas Neumann constant source solution is also presented in 
the figure to illustrate the well-known thin airfoil instability 
associated with source solutions. To account for the lift of the 
supercritical airfoil, the usual approach of adding a uniform 
vortex density around the airfoil was implemented in the Douglas 
Neumann method. It is noteworthy that for this example much of 
the error associated with the source solution could have been 
eliminated by replacing the uniform vortex distribution by a 
distribution parabolic in terms of surface distance (Reference 22) . 
For certain geometries, the parabolic distribution allows weaker 
source strengths and therefore improves the prediction accuracy. 

In order to determine sensitivity to coordinate noise, the 
upper surface velocity distribution was calculated for the non- 
lifting nearly flat plate of Figure 23. Ideally, the smooth upper 
surface should be insulated from the jagged lower surface and the 
calculated results should agree closely with the solution for 
a nonlifting flat plate at -90° incidence. Such is the case 
for all three Green's identity solutions, which are nearly 
identical. Again, the source solution reflects the thin airfoil 
instability . 

The final example explores the limiting behavior of a 
symmetrical airfoil at zero incidence as thickness ratio 
approaches zero (^figure 24) . The perturbation velocity for 
the low order potential boundary condition method and for the 
high order flow tangency method share approximately the same 
average prediction error, which is much less than the error from 
the low order flow tangency solution (figure 25). As a function 
of panel density, the vortex density at 50% chord for all three 
methods approaches the analytical linearized theory value as 
panel density increases (figure 26). However, only the vorticity 
from potential boundary conditions exhibits good accuracy. 
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Figure 23. Noisy Flat Plate at —90° Incidence 

(Cg = 0) 
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Figure 25. Thin Symmetrical Airfoil at Zero Incidence 




Figure 26. Thin Symmetrical Airfoil Vorticity at 50% Chord 



From the above examples , one reaches the heretofore unex- 
pected conclusion that, even in a low order formulation, the 
application of internal potential boundary conditions provides 
prediction accuracy that is in most practical cases equivalent to 
a higher order formulation. This applies only to the velocity 
calculated by equation (20) ; velocity calculated by summing the 
product of velocity influence coefficients and singularity density 
is significantly less accurate. The preliminary explanation is 
that the increment in potential between adjacent control points 
is nearly independent of both source gradient and surface cur- 
vature effects for a wide range of geometric shapes. Therefore, 
the prospect exists for obtaining higher order prediction 
accuracy without the additional complexity of calculating 
higher order corrections to the influence coefficients. 

Accuracy of Concave Corner Solutions 

On typical aircraft configurations, the intersection of 
the wing and fuselage generates a sharp concave corner in the 
cross-flow plane. The existing Douglas Neumann surface source 
method predicts a velocity distribution that increases without 
bounds as the corner is approached, rather than the finite 
velocity limit corresponding to the correct inviscid flow 
tangency solution. This characteristic of the Douglas Neumann 
method is discussed by Hess in reference 23 for the special case 
of two-dimensional flow. 

In order to ascertain whether implementation of Green's 
identity offers the prospect of improving prediction accuracy 
near sharp concave corners, an investigation was conducted for 
the simple two-dimensional geometry of figure 27. This 
geometry was initially used by Hess in his study of source 
method prediction characteristics. He compared the numerical 
and exact analytical solutions near the concave corner at S = 0. 
The present study repeats his procedure with the inclusion of 

predictions by the higher order Green's identity panel method 
using flow tangency boundary conditions. 
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Figure 27. Simple Geometry with a Concave Corner 

Hess provided the solutions from his low order Douglas 
Neumann source method. The geometry of figure 27 modeled by 
Hess incorporates 56 equal length panels between S = 0 and 
0.30 and 66 equal length panels between S = 0.30 and 1.00. 
Identical paneling was used on the lower half of the geometry. 
For the Green's identity panel method, we used essentially the 
same panel distribution, but only one-half the panel density 
in order to avoid increasing the dimension limits of the pro- 
gram. 

Both the symmetric (0° a) and anti-symmetric (90° a) solu- 
tions were investigated. At any intermediate angle of attack, 
the solution is simply a linear combination of these two 
solutions. In the vicinity of a concave corner, Hess states 
that the exact velocity varies as V ^ for symmetric flow and 
V S for anti-symmetric flow. 

The calculated velocity distribution for symmetric flow is 
presented in figure 28 as a function of surface distance S. 

The Douglas Neumann and Green's identity solutions are virtually 
identical except near the concave corner, where the former blows 


50 



up. On the other hand, the Green's identity solution is seen 
to approach stagnation at the corner , as does the exact analyti- 
cal solution. The detailed behavior of the Green's identity 
solution very near the corner is evident in figure 29, where 
velocity has been plotted versus surface distance on a log-log 
scale. Notice that the numerical solution agrees well with the 
exact solution (V S ) except within two panel lengths of the 
corner, where the velocity error is approximately 0.0 01 Voo . 

This error is insignificant and would not be noticeable if a 
logarithmic scale were not used. 



Figure 28. Calculated Velocity Distribution for Symmetric Flow 
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Figure 29. Velocity Distribution Near A Concave Corner 

(Symmetric Flow) 


Similar graphs for anti-symmetric flow are presented in 
figures 30 and 31. Here the discrepancy between the two methods 
is much smaller, although the Green's identity approach still 
provides a significantly better velocity distribution near the 
concave corner. 







Figure 30. Calculated Velocity Distribution for Anti-Symmetric Flow 







The concave corner example was repeated for the low order 
internal potential approach, and the results are similar to the 
solution for the higher order flow tangency Green's identity 
formulation. 

It is concluded that for two-dimensional flow. Green's 
identity provides an improved flow model in the vicinity of a 
concave corner when compared to the Douglas Neumann solution. 

It is expected that this benefit will carry over to the wing- 
fuselage intersection region in the modeling of an aircraft. 

Inverse (Design) Capability 

The approach of reference 17 was used to formulate the 
complete partial derivative distribution of surface velocity 
with respect to surface geometry, for the low order, internal 
potential boundary condition formulation of Green's identity. 
Then, by prescribing an arbitrary change in surface velocity 
distribution, the program determines the corresponding first orde 
change in geometry by solving a system of linear equations. By 
iteration, the program designs the multielement airfoil geometry 
having a specified velocity distribution on one or more elements. 
For each element to be designed in a multi-element airfoil 
system (up to five elements) , the following steps are involved: 

(1) The user prescribes a design pressure or velocity 
distribution around the surfaces of one or more of the 
various airfoil elements. The geometry of the remain- 
ing elements will be fixed. 

(2) The user prescribes a starting geometry to initialize 
the calculations and the location of one point on each 
element to be fixed in space, such as the trailing 
edge . 

(3) The program solves the direct problem for the starting 
geometry in order to determine the change in velocity 
distribution required to achieve the prescribed values. 

(4) The program calculates the rate of change of surface 
velocity with respect to an arbitrary change in surface 
angle distribution. The element perimeter remains 
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(5) 


fixed. If the tangential component of velocity at the 

t h 

control point of the i^ panel is designated Vt-, and 

th ^ 

if the surface angle of the j panel is designated 0j , 
then the array is calculated where 



(23) 


The change of surface angle distribution is calculated 
in accordance with the prescribed velocities and the 
following first order expression. 


AV^ = I (A^. A0.) (24) 

i j J 

(6) The geometry is corrected by the program, and steps 

(3) - (5) are repeated as a series of iteration cycles. 

The most difficult and important step in formulating the inverse 
capability is to generate the matrix A^^j . It is noted that all 
terms were incorporated in deriving the partial derivative, 
including singularity strength changes and the displacement of 
panels j+1, j+2, etc. corresponding to the surface angle change 
A0j. The corresponding singularity strength changes are obtained 
by a first order expansion to the boundary condition equation. 

A typical inverse solution requires five iteration cycles, 
where each cycle requires approximately four times the computing 
time of a direct problem solution. On the CBC CYBER 173, a 
typical two-element airfoil inverse problem with seventy panels 
uses 20 seconds computing time per cycle. This compares with 
150 seconds per cycle for the earlier least squares method 
(reference 17) . 

Two examples are discussed herein to illustrate inverse 
problem solution convergence characteristics. 
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The objective of the first example is to design a circular 
cylinder by using a nearly flat plate for the starting geometry 
(Figure 32). The exact analytical surface velocity distribution was 
prescribed, and the converged solution geometry of figure 32 
was obtained after four iteration cycles. The panel endpoints 
are within a maximum distance of 0.002 x radius of lying on a 
circle. The complete partial derivatives of velocity with 
respect to surface angle change are necessary but not sufficient 
for obtaining convergence about the periphery of this example. 

The use of mild combined source-vortex singularities is also a 
factor. To illustrate, the example was repeated, but this time 
only vortex singularities were used to induce the flow field, in 
accordance with the analysis method of Woodward-Dvorak (reference 
11). The geometry never converged (figure 15) but oscillated 
+30° in the leading edge region from one iteration cycle to the 
next . 

The second example demonstrates the inverse solution for the 
two-element Williams airfoil presented earlier in figure 16. The 
simple starting geometry shown in figure 33 was used to initialize 
the calculations, and the exact surface velocity distribution of 
Williams was prescribed on both elements. The geometry converged 
and agreed with the target geometry to within a tolerance of one- 
tenth of one degree surface angle on all panels in five iterations. 

This approach to the design problem is suitable for exten- 
sion to 3-D. The 2-D inverse method based on the Green's 
identity formulation has demonstrated low cost combined with 
consistent accuracy and stability. The low calculation cost will 
be even more important in a 3-D method. The effort of developing 
a 3-D method is reduced since the modified 3-D analysis program 
uses the same low-order singularities and internal potential 
bomdary conditions as the 2-D inverse method. The modified 3-D 
program can be used directly for the analysis steps. The inverse 
capability would play an obvious role in future high-lift device 
design, and would also be valuable in analyzing flow fields with 
strong viscous-inviscid interactions. 
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Selected Formulation 

Overall, the research into two-dimensional panel method 
formulations indicates that only through the application of the 
combined source-doublet distribution of Green's identity will 
one consistently obtain numerically stable calculations in both 
analysis and design modes without any significant increase in 
computational expense compared to either a source or doublet 
only formulation. With higher order corrections, flow tangency 
boundary conditions consistently provide accurate predictions. 
Competitive accuracy can be obtained for a wide range of shapes 
without such corrections if internal potential boundary condi- 
tions are applied. 

For initial development of a 3-D Douglas Neumann Program 
modification, the internal potential boundary condition formula- 
tion has the advantage of using the low-order source panel 
singularities already in the program, while the flow tangency 
boundary condition approach would require the additional effort 
of incorporating the high-order, curved panel singularities 
still under development by Douglas Aircraft Company. The internal 
potential method is also a consistent analysis method for extending 
the 2-D inverse method to a 3-D capability. Therefore, the inter- 
nal potential formulation was selected for the 3-D program modifi- 
cation . 
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NUMERICAL SOLUTION FORMULATION 


The formulation is presented as a modification to the exist- 
ing Douglas Neiimann Program (reference 2) to improve the effi- 
ciency and reliability of predicting component interference and 
high lift characteristics of wing-body configurations. The pri- 
mary modification is the substitution of the mild combined source- 
doublet distribution of Green's identity for the fundamental 
source-only distribution of the existing program. The advantages 
associated with combined source-doublets were explained in terms 
of classical theory in the section on Potential Flow Theory and 
were demonstrated for two-dimensional flow in the examples of the 
section Research on Green's Identity Formulation. 

The objective is to couple contemporary numerical techniques 
with the classical theory in order to provide a reliable predic- 
tion tool. It is assumed that one or more bodies are immersed in 
a steady, incompressible, inviscid stream of velocity V„ . The 
bodies can be either lifting or nonlifting, v/hich means that wakes 
can be either included or omitted. Consistent with the existing 
program, the present formulation assumes solid body tangency 
boundary conditions. In accordance with the classical theory, it 
would be a simple matter to extend the formulation to allow either 
arbitrary normal velocity boundary conditions or Dirichlet pre- 
scribed potential distribution. The aim is to calculate flow pro- 
perties at the solid body surfaces, particularly the pressure, 
which is integrated to provide force and moment. 

The selection of the modeling is dictated by the following 
guidelines : 

1. If the body geometry approaches the shape of an infinite 
aspect ratio wing with uniform cross-section, then the 
numerical doublet model should reflect the piecewise 
linear vortex representation employed by the two-dimen- 
sional Green's identity panel method of the preceding 
section . 

Experience in solving two-dimensional problems demon- 
strates that such a vortex model provides good accuracy 
and numerical stability, and it would be unfortunate to 
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sacrifice such characteristics in the three-dimensional 
formulation. 

The doiiblet distribution model should be as nearly con- 
tinuous as possible. Discontinuitites in doublet density 
correspond to concentrated vortex filaments which coun- 
teract the mild velocity gradients associated with 
distributed surface singularity methods. 

3. The entire doublet distribution should be described by 
a set of scalar parameters with approximately the same 
number of members as there are panels. These parameters 
are the basic unknowns to be determined. Determination 
of additional parameters could lead to excessively large 
computational expense, since computing time varies as 
the third power of the number of unknowns . 

There are four basic steps in the solution process, (1) geo- 
metry definition, (2) influence function calculation, (3) solution 
to linear simultaneous boundary condition equations for the un- 
known singularity density distribution, and (4) calculation of 
surface flow properties. These steps correspond to a direct 
application of the theory, except that the continuous integral 
theoretical relationships are discretized to allow numerical 
solutions to arbitrarily complex problems. Prediction error is 
defined as the difference between the exact analytical solution 
and the numerical calculations. Because no small disturbance 
assumptions are employed, prediction errors result from the 
discretization process only. Therefore, one could generate 
arbitrarily accurate numerical predictions at the expense of the 
increased computing time commensurate with the introduction of 
additional unknowns. 

The geometry definition involves replacing the actual contin- 
uous boundaries by a set of trapezoidal panels. In the existing 
program of Ref 2, each panel is identified by one of four categories 
(1) nonlifting, (2) lifting, (3) internal lifting, and (4) wake. 
Nonlifting source only panels are typically distributed on 
fuselage- type surfaces while the lifting panels model wings and 
aerodynamic control surfaces. The internal lifting panels are 



designed to pass through fuselages and connect the exposed wing 
roots in order to provide fuselage lift-carryover. In the present 
formulation, the combined source-doublet distribution panels 
allow a natural carryover of flow properties from wing to 
fuselage surfaces. Therefore, the first three panel catego- 
ries of the existing program are replaced by the single category 
"body panels" . In the same sense that the theory recognizes no 
fundamental distinction between fuselages, wings, and tails, the 
body panels are used to represent any solid flow boundary. As 
in the existing program, wake panels are used to model lifting 
configurations . 

The continuous surface source and doublet distributions of 
the theory are respectively modeled by piecewise constant and 
piecewise quadratic distributions on the trapezoidal panels. 

Seven parameters describe the complete singularity density 
distribution on each panel. One parameter is the uniform source 
density. The remaining six are the coefficients of the six terms 
in the expression for quadratic doublet density. The potential 
and velocity induced at any field point by the panel singularity 
distribution can be expressed as the product of each of the 
seven parameters with its appropriate influence function. Each 
influence function depends solely on panel geometry and field 
point location and is written in closed analytical form. One 
influence function describes the effect of a unit strength uniform 
doublet distribution, another describes the effect of a unit linear 
doublet distribution, etc. 

Of the six parameters that describe the doublet distribution 
on a panel, only one can be considered an independent variable 
reserved for boundary condition satisfaction. The remaining 
five are fall-outs whose values are dictated by a least squares 
surface fit through neighboring control points. 

There is approximately one boundary condition control point 
allocated per panel. Because Green's identity is employed, the 
source density for each panel is a function of local geometry 
only (eq. 16) . It remains to determine the doublet density. 

Rather than prescribing Neumann boundary conditions directly. 
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the theoretically equivalent zero internal perturbation potential 
boundary conditions are applied. In two-dimensional flow, such 
potential boxmdary conditions consistently produce very good 
prediction accuracy even though panel curvature and source gra- 
dient effects are neglected. Furthermore, the two-dimensional 
design boundary conditions for inverse problems were easy to 
derive when coupled with internal potential conditions and led to 
exceptionally rapid, stable convergence. Based on the experience 
that three-dimensional numerical characteristics should reflect 
those of two-dimensional flow, the selection of internal potential 
boundary conditions was made. 

With the aid of the influence functions and doublet con- 
tinuity considerations at panel edges, the internal potential 
boundary conditions are established as a system of linear equa- 
tions with the doublet strength parameters as unknowns . There 
is essentially one unknown per panel. Solution to the system 
renders the complete singularity distribution known, and it is 
then a simple m.atter to calculate the net induced velocity or 
potential at any field point. Pressure coefficients at body con- 
trol points are determined from Bernoulli's equation (6), and 
the resultant force and moment is determined by integration 
under the assumption that the pressure on each panel is uniform. 

The formulation is presented in detail below. 

Geometry Panel Modeling 

All boundary surfaces are divided into continuous regions 
designated sections, as illustrated in figure 34. Each solid 
body in the flow can be described by one or more sections, and 
the division is at the prerogative of the user. For example, a 
section could be the upper surface of a wing, a fuselage, or 
both the fuselage and wing simultaneously. The only restriction 
is that a single section must contain either body regions or wake 
regions exclusively and not both simultaneously. This deviates 
from the original program but provides a more consistent formu- 
lation. It is assumed that all body slopes within a section are 
continuous; therefore, section edges should be aligned with 
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any slope discontinuities such as wing trailing edges and wing- 
fuselage junctures. 



Each section is subdivided into panels in the same manner 
employed by the existing Douglas Neumann Program. The section 
is described by a total of NxM points, where each point is iden- 
tified by the pair of indices i and j (1-ifN, lijfM). 

See figure 35. The points describe a set of (N - 1) x (M - 1) 
panels, where each panel is identified by the pair of indices 
ip and jp ( 1 £ ip < N-1; 1 < jp < M-1) . See figure 36. Panel 
(ip,jp) is defined by the four corner points (ip, jp) , (ip, jp+1) , 
(ip+1, jp+1) , and (ip+1, jp) . The order of sequencing is such 
that the positive normal direction (pointing into the fluid) is 
out of the plane of the paper. In other words, the positive 
normal direction is in the sense of the increasing i-direction 
crossed with the increasing j-direction. Any inertial (x,y,z) 
Cartesian coordinate system is satisfactory for describing the 
geometry. The free stream velocity is assumed to have unit 
magnitude and to be described by the following components: 
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( 25 ) 


where e 


X 




V = a^e + a e,, + a e^ 
CO XX y y z z 


e are unit vectors in the x, 
z 


and z - directions. 


respectively . 

The four corner points of a panel in general describe a 
nonplanar quadrilateral. For compatibility with existing 
influence functions, adjustments to the corner points 
are made in order to generate the trapezoid that is most nearly- 
described by the original four points. The procedure is iden- 
tical to that of the existing program (reference 2) . 

Consider the arbitrary panel of figure 37. For clarity, the 
four corner points are identified by indices k = 1, 2, 3, 4. 

The adjustments to determine the trapezoidal panel involve 
making the line between points 1 and 2 parallel to the line be- 
tween points 3 and 4. If the unadjusted coordinates of corner 
point k are (^j,/ yj^ ^ ^k^ ' then define vector as 


E X, e -t y, e + z, e„ 
k X -'k y k z 


(1 - k - 4) 


(26) 



Figure 37. Adjustment of the Input Points to Form a Plane 
Trapezoidal Panel 
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and, using subscripts "F" 
define 


and "S" to designate first and second. 



)>(27) 
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-y -y 

The weighted average of and Pg is used to define the 
direction of the parallel sides of the trapezoid, which is also 
selected as the C~d.irection of the trapezoidal panel coordinate 
system 


-y 




(28) 


If the adjusted corner points are identified by an asterisk *, 
then the following definitions are applied: 
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It is noteworthy that the midpoints and lengths of line segments 
1-2 and 3-4 remain unchanged after the adjustment. 

The normal direction is defined as the direction of the 
?-axis of the panel coordinate system 


^ (?^* - ?2*) X (?3* - ?^*) 

0 _ __ _ _ _ 

^ |(P^* - P2*) X (Pg* - P^*) I 


(30) 
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The third Cartesian coordinate of the panel system is designated 
n and is assigned the unit vector e^ where 


A 3 X 3 transformation matrix a^^ is generated from the above 
definitions such that 


{e 






e^}= [ai_] 



(32) 


The centroid of the trapezoid is taken as the origin of the panel 
coordinate system ii, n / ?)• See figure 38. 


(ip, jp+1) 



Figure 38. Panel Coordinate System 


Several important geometric parameters associated with each 
trapezoidal panel such as area, maximum diagonal, etc., are 
calculated for future reference. These are described in detail 
in reference 2 and are not repeated here. 
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Panel Singularity Distribution 

The theoretical singularity distribution on the actual bound- 
ary surface is modeled by a uniform source density and quadrat- 
ically varying doublet density on each trapezoidal panel. As 
part of this contractual effort, John Hess of Douglas Aircraft 
Company furnished a computer program code to calculate both the 
velocity and potential induced at any specified field point by 
a uniform source density and arbitrary quadratic doublet density 
distribution on an arbitrary trapezoidal panel. In order to 
generate reasonable doublet continuity modeling at panel edges, 
an approach developed by Boeing (reference 16) has been applied. 
This approach establishes linear relationships between the coef- 
ficients of the quadratic doublet distribution on each panel with 
the coefficients on the adjacent eight panels. An alternate 
approach which is believed to provide better continuity proper- 
ties at the expense of slightly increased computing effort has 
been formulated as part of the contractual study, but has not yet 
been coded for a computer. These developments are discussed 
below. 

Influence Function Formulas - Consider a trapezoidal panel 
with the geometry of figure 38. The uniform source density is 
a and the quadratic doublet distribution p(C,n) on the panel is 
described as follows: 

where the coefficients Mqq/ •••/ Vq 2 arbitrary. The formulas 

furnished by Hess provide the induced potential and velocity at 
an arbitrary field point n^, • In Hess' formulas, the 

coordinates are designated (x, y, z) . It is impor- 

tant to interpret the coordinates as being in the panel system ^ 

(5, n ? c) • The term (j)g is defined as the potential induced by 
a uniform source density of unit strength, while refers to 

the potential induced by the doublet density distribution 
y(C,n) = n^. Then the potential ^ induced by all the singu- 

larties on the panel is expressed as follows: 
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The velocity is generated by taking the gradient of <t> with 
respect to ( 5 ^, ' i.e., with respect to (x, y, z) in Hess' 


formulas . 
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Analogous expressions exist for |^— and . 

The required influence functions are , *^00' •'*' *^0 2' 

3<f>g ^*^00 ^*^02 . Each of these twenty-eight functions 

^x ' "^jT • ' 3z 

depend solely on the coordinates of the field point and the four 
corner point coordinates of the trapezoidal panel , 

Some of the functions are provided 
in reference (2), along with intermediary functions of the four 
corner points. Using the same definitions from pages 77-84 of 
reference (2), the remaining influence functions are presented 
below. It is important to realize that Hess ' definitions of 
potential and velocity influence functions differ from those of 
this report by the factors -4 tt and +4ir , respectively. To avoid 
confusion, Hess' formulas will be designated by an asterisk, 
where 

(p * B — 4 IT (f) 


V* B +4ttV 
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convert 
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from 

the 

panel 

coordinate system to 

the body system 

of figure 

34, the transformation 

matrix 

is employed. 










It 

is to 

be 

noted that 

the 

influence 

functions are exact for 


any field point location. Approximate intermediate-field and far- 
field formulas were not furnished by Hess and have not yet been 
incorporated. The advantage of such approximations would be 
a significantly reduced computing effort at no appreciable loss 
of accuracy for points not in the immediate vicinity of the 
panel. For example, in the far field representation the effect 
of the uniform source density a on the trapezoidal panel of area 
A is represented by a point source of strength aA positioned at 
the panel centroid. The formulas for a point source are, of 
coiirse, far simpler than those for a distributed surface singu- 
larity. 
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Doublet Distribution Surface Fit Approaches - On each panel 
the doublet distribution is described by the six quadratic 
coefficients viqq/ • • • » ^02' Several approaches were consi- 

dered for matching the distributions on neighboring panels in 
order to accurately model the analytic doublet distribution and 
to minimize doublet discontinuities at panel edges. An approach 
developed by Boeing (Reference 16) was selected for implementa- 
tion in the Douglas Neumann Program modifications. The essence 
of this approach involves passing a least squares quadratic 
through the boundary condition control point of a panel and 
through the control points of the adjacent eight panels. An 
alternate Boeing approach for generating improved doublet con- 
tinuity characteristics was also coded and tested, but proved to 
be numerically unstable if the panel network is not composed of 
nearly straight generator lines. A new approach that improves 
both the continuity and accuracy of the quadratic doublet 
representations without introducing numerical instabilities 
was developed under this study but has not yet been coded for 
the computer. A discussion on the various formulations is 
presented below. 

The numerical requirements for achieving continuity can be 
illustrated by considering a broken line segment model of a 
continuous function y = f (x) . See figure 39. Over each inter- 
val i, the model is described by y = a^x + b^ , where the 
coefficients a^ and b^ are analogous to the six quadratic 
coefficients VJ qq f • • • ' ^ 02 ^ panel in the surface doublet 

distribution. Consistent with the continuity of function f(x) , 
the broken line segment model is required to be continuous at 
interval endpoints. If x^ is the endpoint between intervals 
i and i+1, continuity is imposed by the condition: 


a.x. +b. =a.,,x. +b.,. 
IX 1 x+1 X x+1 


which can be re-expressed as: 


a . , T = a . - 
x+1 X 


(bi^l - b,) 


X . 
X 


( 42 ) 
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Figure 39. Broken Line Segment Approximation to Function f(x) 


Regardless of the function f(x), equation (42) expresses the 
continuity constraint that is a linear function solely of 

a^, b^, and other words, there is approximately one 

free parameter per interval for adjusting the broken line seg- 
ment model to most nearly match the curve f (x) . The remaining 
coefficients are entirely dictated by continuity. 

The necessity for minimizing doublet discontinuities is 
apparent upon examination of a simple numerical example. Con- 
sider a two-dimensional doublet distribution described by the 
equation 

3 

y = S 

where y is doublet density and s is surface distance. A cubic 
has been selected because it cannot be modelled exactly by 
quadratic curves- Now suppose piecewise quadratic interpolating 
curve fits are made to the doublet distribution between integral 
values of s. That is, the endpoints are at s = ..., -2, -1, 0, +1, 
+2, ... Of the many possible types of quadratic fit, two are 

illustrated in figure 40 over the range 0 < s < 2. The first 
curve fit (A) is continuous with continuous slopes at panel 
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endpoints, equivalent to the piecewise linear vortex distribu- 
tion employed in the method of references (19-21). in this case, 
the quadratic fit was determined such that the exact slope is 
attained at all panel endpoints and such that the curve fit 
passes through the exact value of doublet density at s = 1 . The 
second quadratic curve fit (B) was established at each interval 
by passing a parabola through the exact doublet density at the 
interval midpoint and at the midpoints of the two adjacent inter- 
vals. At this stage it is difficult to predict which type curve 
fit would be most accurate for use in a numerical surface 
singularity method. Whereas the discontinuity in curve fit B at 
s = 1 is obviously nonexistent in the exact distribution, such 
discontinuities do allow quadratic fit B to attain a more accurate 
average doublet density in each interval than curve fit A. The 
analytic gradient of the doublet distributions provides the 
equivalent vortex density curve fits (figure 41) . For curve fit 
B, the doublet discontinuity at s = 1 is reflected by a delta 
function in the vortex distribution. 



S - Surface Distance 


Figure 40. Two Curve Fits to a Doublet Distribution 
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Figure 41. Vortex Distribution - Corresponds to Figure 40 

Neither the doublet nor vortex distributions reveal the sig- 
nificant numerical consequences of the type of curve fit. How- 
ever, the corresponding normal velocity distribution is highly 
informative. For both curve fits and for the exact doublet 
distribution., the normal velocity induced by the section of dou- 
blet distribution between s = 0 and s = 2 was calculated analyt- 
ically. A finite vortex filament was included in each case at 
s = 2 to counteract the doublet discontinuity magnitude v = 8. 

Plotted as a function of surface distance s in figure 42, 
the normal velocity distribution corresponding to curve fit B 
is seen to be the more accurate in the vicinity of s = 1/2 and 
s = 3/2, which would correspond to control point locations in a 
panel method. However, the effect of the concentrated vortex 
at s = 1 leads to numerical instabilities for curve fit B. The 
consequences are quite clear. No matter what type panel spacing 
is used or how close a control point is to a panel edge, curve 
fit A will provide a reasonable approximation to the actual exact 
normal velocity curve. On the other hand, curve fit B is obvi- 
ously in error unless the control point is very close to the 
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panel midpoint. Therefore/ erroneous calculations could result 
if prescribed normal velocity boundary conditions were used with 
curve fit B. Although it is less obvious, the application of 
internal potential boundary conditions would not eliminate the 
problem. In fact, for a thin wing with upper and lower surface 
panels sharing a common chord plane projection, the use of in- 
ternal pobential boundary conditions approaches equivalency to 
lifting surface theory as wing thickness vanishes. It is con- 
cluded that minimizing doublet discontinuity is highly beneficial 
to the modelling of a doublet distribution. 



Figure 42. Normal Velocity Distribution - Corresponds to Figure 40 

If a three-dimensional doublet representation is to satisfy 
the three guidelines presented in the introduction to this 
section, it is necessary that the doublet density on each panel 
vary as a second order polynomial in terms of coordinates in the 
plane of- the panel. The objective is to adjust the coefficients 
such that boundary conditions are satisfied at approximately one 
control point per panel with as smooth a doublet distribution 
on the panel network as possible. 
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A schematic of the selected control point locations is 
furnished in Figure 43. For each trapezoidal panel there is 
one control point at the centroid. Additional control points 
are placed along the edges of a section at the trapezoidal panel 
side midpoints and at section corners. The purpose of the 
additional points is to generate doublet continuity at adjacent 
section edges without requiring two-sided surface fits to the 
doublet distribution across the edges. 



Figure 43. Schematic of Panel Control Point Locations on a Section 


The entire doublet distribution on a section is uniquely 
determined by the set of doublet densities at the control point 
locations. On each panel the distribution which has been 
selected for use in this study is generated by a least squares 
quadratic fit through nine neighboring control points. For 
example, consider any panel (ip/ ^p^ ' which is schematically 
illustrated in Figure 44. The doublet density at the panel 
centroid control point and at the eight adjacent control points 
is identified by subscript k (1 ^ k £ 9) . If panel (i^f j„) is 
on the edge of a section, then some of the adjacent control 
points will be on panel edges. The following table provides 
the conversion between control point index k of Figure 44 and 
indices (i, j) of Figure 43. 
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Figure 44. Indexing Convention for m in Vicinity of Panel (ip, jp) 


TABLE II - CONTROL POINT INDEX (i, j) 
[Panel Index (i , j )] 

ir 




1 


^P 


^P 

2 


^P 

^P 

+ 1 

3 


i 

P 

^P 

+ 2 

4 

i 

P 

+ 1 


^p 

5 

i 

P 

+ 1 

^P 

+ 1 

6 

^P 

+ 1 


+ 2 

7 

"p 

+ 2 


^P 

8 

i 

P 

+ 2 

^P 

+ 1 

9 

i 

D 

+ 2 


+ 2 


The second order polynomial u(C,n) of Equation (33) is the 
doublet distribution on panel (ip/ jp) • The six coefficients 
Uqq,. . . , Pq 2 determined by minimizing the following error 

function E 
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E 


2 


(43) 


9 

= I 
k=l 


Pj,]} 


(C]^» rij,) are the coordinates n) of the kth control point. 

For control points k that are either at the centroid or on the 
edge of panel (i_, ^ very large weight factor W,>>! is 

selected. Otherwise, is chosen as unity. This weighting 
matches the function p) to the values at control points 

on the panel and provides approximate matching at the remaining 
control points. 

Equation (43) is minimized with respect to each coefficient 
Uqq,. . . , Vq 2 * This leads to a system of six linear equations. 

If the six coefficients y^^, y^^, y^^, y^^, y^^ 

are respectively identified by 62 /- . . , gg, then the 

solution to the system of equations can be expressed in the 
following form: 

9 

^ ^k (ll<tl 6 ) (44) 

k=l 

The array is determined by matrix inversion and depends 

solely on the panel geometry, not on the values yj^. Therefore 
array can be determined before solving for the doublet 

strength . 

The above surface fitting approach for establishing the 
doublet distribution is equivalent to method "B" in the two- 
dimensional example of Figures 40-42. Admittedly, the resultant 
doublet continuity properties at panel edges are not ideal. An 
alternate approach developed by Boeing in an unpublished report 
reduces to method "A" for the special case of two-dimensional 
flow. This alternate approach was coded under the present study 
and is available in the modified Douglas Neiamann Program. For 
cases in which the paneled geometry of a section is defined by 
(nearly) straight generator lines, the alternate approach is 
accurate and has good continuity properties. Otherwise, however, 
the calculated array analogous to B , tends to be numerically 
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unstable. Therefore, the former approach was selected and was 
used in all the three-dimensional examples of this report. 

A third approach, previously investigated by Boeing, involves 
reducing the discontinuities at panel edges after the doublet 
distribution has been determined. The computed doublet densities 
at panel edges from the initial doublet distribution are averaged 
to generate a new, nearly continuous distribution. Then a least 
squares second order fit is passed through the new edge values. 

The drawback is that the discontinuities should be eliminated 
prior to satisfying boundary conditions, not after. The signifi- 
cant numerical instability associated with doublet discontinuities 
is that boundary conditions become unrealistically sensitive to 
panel geometry and control point location. 

An improvement to the third approach has been developed 
under this study, but has not been coded at this time. The 
essense of the present improved approach involves minimizing doub- 
let discontinuities at eight peripheral points on each panel prior 
to solving for the unknowns. This is possible because the doublet 
distribution on each panel is a linear function of the unknowns, 
where the unknowns are the doublet density at the adjacent con- 
trol points of Figure 43. Four of the eight peripheral points on 
a panel are the corners and four are the side midpoints (Figure 45) 
In order to assure that adjacent panels in a section will share 
identical peripheral point locations at a corner or edge, the 
eight point locations of Figure 45 are to be interpreted as lying 
on the sides of the non-planar quadrilateral panel which was de- 
fined prior to the trapezoidal approximation. Doublet discontin- 
uities at panel edges are minimized by attempting to match the 
doublet distributions of adjacent panels to a common desired value 
at each of the peripheral points shared by the adjacent panels. 
Details of the present improved approach are discussed below. 

A second order quadratic fit y(5, n) is to be established 
on each panel 

y( 5 ,n) = 6^ + 62C + 33 n + 3 ^ 5 ^ + n + 3g ,45, 
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Figure 45. Panel Edge Locations for Doublet Matching 

For each panel, the coefficients (1<_S.^6) are determined by an 
exact matching of y(5,n) to the unknown centroid value and a 
least square error matching to desired doublet values at the 
eight panel peripheral points. For any peripheral points that 
are on the edge of the section, the square error weighting is 
selected to make the matching exact. If is the unknown 
centroid doublet density of the panel and ypj^ (l£k_^8) identifies 
the desired doublet density at the eight peripheral points, then 
the method of least squares generates a relationship in the 


following form: 



The coefficients and are functions of only the panel 

corner coordinates (prior to the trapezoidal approximation) . 

By expressing ypj^ (l£k£8) as a linear combination of the 
unknown doiablet densities at several neighboring control points 
and then substituting the expression into Equation (46), the 
coefficients 3 ^ of a panel will be in the following form; 


= I B. y 
£m m 
m 


(1<A<6) 


where y refers to the unknown doublet densities at control 
points in the vicinity of the panel. 
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To determine as a linear combination of the unknowns, 

an interpolating least squares surface fit is passed through the 
control point doublet densities in the neighborhood of the panel 
peripheral point k. For best interpolation accuracy, the 
greatest least squares weighting should be assigned to control 
points closest to point k. 

The array of Equation (47) is a function only of the 

corner point coordinates of neighboring panels and is determined 
prior to solving for the unknowns The limits of m in 

Equation (47) depend upon the number of neighboring control 
points used in the interpolation for ppj^- The limits l£mj£2 5 are 
expected to be adequate. 



Figure 46. Three Quadratic Curve Fits to a Doublet Distribution 

The present approach is expected to improve both the 
accuracy of the doublet distribution and the continuity proper- 
ties, at the expense of slightly greater computational effort. 
However, it is expected that the increased effort will be 
insignificant compared to what is required in establishing 
influence coefficients and solving the system of linear boundary 
condition equations. 
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The numerical characteristics are aptly illustrated by 
re-examining the simple two-dimensional example of Figures 40-42 
with the inclusion of the present approach (curve C) . It is 
noted that in terms of doublet density (Figure 46) , vortex 
density or doublet gradient (Figure 47) , and induced normal 
velocity (Figure 48) , the present method most nearly matches 
the analytic curve. The significance is revealed in the normal 
velocity distribution of Figure 48. Regardless of where one 
might select a boundary condition control point, the boundary 
value is sufficiently close to exact to suppress unwarranted 
numerical instabilities. Therefore, panel geometry and control 
point location can be expected to have only a minor effect on 
computed results, which is compatible with simplified user 
requirements . 



Figure 47. Vortex Distribution - Corresponds to Figure 46 
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Figure 48. Normal Velocity Distribution - Corresponds to Figure 46 

Boundary Conditions 

The indirect approach of applying internal potential boundary 
conditions is applied in view of the good accuracy and stability 
exhibited in the two-dimensional examples presented earlier. 

The two steps to satisfaction of the boundary condition are to 
set the source density on each panel equal to the negative of 
the local free stream normal velocity component and subsequently 
to determine the doublet distribution corresponding to zero 
internal perturbation potential at selected control point loca- 
tions. As explained in the section Potential Flow Theory, such 
an approach is theoretically equivalent to prescribing solid 
body flow tangency conditions. 

Several other approaches for applying boundary conditions 
are conceivable. For example, the net velocity flux through a 
control surface bounded by four neighboring panel centroids 
could be prescribed as zero. This would require analytical 
integration to formulate velocity flux influence coefficients 
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for three-dimensional geometries, but would reduce to the suc-r 
cessful stream function approach of Oellers (Ref 10) for two- 
dimensional flow. However, the present internal potential 
boundary condition approach is believed to provide the best 
compromise with respect to accuracy, numerical stability, compu- 
tational efficiency, and availability of influence functions. 

The source density on each body panel is calculated from the 
dot product of the free stream velocity with the local normal 
direction e^, in accordance with equation (16) . 


a 


V 


^“x^31 “y^32 “z^33^ 


(48) 


For wake panels the source strength is prescribed to be zero, 
corresponding to the assumption of negligible viscous displace- 
ment effects. 

A schematic of the control point locations was presented 
earlier in Figure 43. The edge control points do not lie pre- 
cisely on the trapezoidal panel edges, but are instead moved a 
few percent of local panel dimensions toward the centroid. 

This repositioning prevents contact with the concentrated 
vortex filament that bounds the edges of a doublet sheet. 

At each control point on a section of body panels, the 
boundary condition is imposed that the internal perturbation 
potential induced by the simultaneous action of all singularities 
is zero. The internal potential is evaluated at c = 0 , i.e., 
on the non-fluid surface of the panel. At the control points 
along one edge of each section of wake panels (Figure 49) , the 
imposed boundary condition is that the total velocity component 
normal to the panel is zero. By aligning this section edge with 
the trailing edge of a lifting body, the Kutta condition of wake 
tangency is satisfied. Zero doublet gradient is prescribed along 
the j -direction of the wake panels, corresponding to the absence 
of wake loading. This eliminates the need to prescribe boundary 
conditions at the other wake control points. 
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Figure 49. Schematic of Wake Boundary Condition Control Points 

The present approach of prescribing zero internal perturba- 
tion potential at edge control points provides for a continuous 
carryover of doublet density from one section to another. On 
each side of the intersection between two sections, the edge 
control point will have the same value of internal perturbation 
potential, zero. Because the difference in potential across the 
intersection is dominated by the magnitude of the local doublet 
discontinuity, the absence of any difference assures doublet 
continuity. The discontinuity in doublet gradient (vortex 
density) is dictated by the change in potential between the edge 
and centroid control points of a panel. In fact, the present 
approach is equivalent to the formulation for predicting the 
vortex discontinuity at sharp corners which was described in the 
section "Research on Green's Identity Formulation." Therefore, 
the present method can be applied to sharp concave corners at 
wing-fuselage intersections. Whereas the absence of a sharp 
corner does not negate the applicability of the approach, it 
would probably be preferable to eliminate edge control points 
in cases in which the body is smooth at the section edge. This 
would require the introduction of special indexing to carry the 


88 



doxiblet distribution surface fit from one section to another. 

It is anticipated that such indexing will be established in the 
future . 

The prescription of zero internal perturbation potential at 
edge control points is believed to be less sensitive to gaps 
and control point location than the Boeing approach of prescrib- 
ing zero normal velocity at the same points (Reference 16) . For 
example, consider a simple two-dimensional problem in which two 
adjacent parallel panels exhibit a slight gap (Figure 50) . 

Assume that there is a unit strength doublet density on both 
panels. The corresponding potential and normal velocity compo- 
nent induced by the two doublet panels is plotted as a function 
of position in Figure 50. Whereas the normal velocity distribu- 
tion is a singular function of both position and gap at the panel 
edges, the distribution of potential is well-behaved on both 
panels. Similar behavior is evident regardless of the angle 
between the two panels. Hence the present approach should be 
consistently reliable even without the introduction of special 
additional treatment at section edges. 



PANEL A 
LENGTH 1.00 

M= 1 



LENGTH 1.00 

M= 1 


Figure 50. Potential and Normal Velocity in Vicinity of Adjacent Panel Edges 
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By combining the influence functions of equations (36)- (41) 
with equations (34) , it is possible to establish a system of 
linear equations relating the unknown set of control point 
doublet densities to the imposed boundary conditions. The number 
of equations and unknowns are equal, and the solution renders 
the complete singularity distribution for all panels known. 

In order to simplify the solution for more than one free 
stream vector V^, three right-hand-sides to the system of 
equations are solved simultaneously, with each right-hand-side 
denoted by a different subscript. 

(49) 

Then for any vector ^ , the solution singularity strengths can 
be generated as a simple linear combination of the solutions for 
^ 002 ^ and 



Calculation of Flow Properties 
At each panel centroid the equality between tangential per- 
turbation velocity component and doublet gradient is employed to 
calculate the total local flow velocity V. In the panel coordi- 
nate system (5 /Ti,?)» 


^)e 
H ^ 


+ «2>®{ + + ®3>®n 


The components in system (x, y, z) are obtained with the aid of 
the rotation matrix [a^j]. 

Pressure coefficient is calculated by Bernoulli's equation. 



00 
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Force and moment integration is performed under the assump- 
tion that the centroid pressure of each panel acts on the entire 
flat trapezoidal geometry. The resulting integration accuracy 
is consistent with the basic solution formulation. 


91 



CALCULATED RESULTS 


Surface pressure distributions were calculated for several 
three-dimensional geometries in order to assess the relative 
merits of the modified program developed under this study and 
the existing Douglas Neumann program. Hereinafter, the modi- 
fied and existing programs are respectively designated "present 
method" and "Hess program." The authors are grateful to James 
Thomas of NASA, Langley Research Center, for his interest and 
effort in the selection and testing of the examples. 

The geometry selection included shapes for v;hich prediction 
accuracy by the Hess program is characteristically good as well 
as shapes for which the program tends to be unreliable. The 
former category includes isolated solid bodies of revolution, 
wings of conventional section geometry, and typical wing-fuselage 
combinations. The latter category includes wings with thin, 
highly loaded trailing edge regions and internal duct flow. 

The first example is intended to reveal whether calculated 
pressure distributions for wings of high aspect ratio tend to 
reflect the 2-D numerical characteristics discussed earlier in 
the section "Research on Green's Identity Formulation." The 
supercritical geometry of Figure 51 was panelled as an unswept 
wing of rectangular planform with constant cross-section and 
aspect ratio 100. Ten equally spaced spanwise strips of panels 
were employed. Two chordwise panel spacing distributions were 
examined, each having a total of 40 panels per spanwise strip. 
Both the spacings , cosine and cos/sinh, are dense in the trailing 
edge region, with the former being the denser of the two. The 
calculated chordwise pressure distributions near the root are 
presented in Figure 51 for the two programs and for the virtually 
exact 2-D conformal mapping solution of Catheral-Sells (Reference 
5) . The present method solution is nearly independent of the 
panel spacing, whereas the Hess program solution is highly un- 
stable. It is important to consider that in the Hess program 
a uniform chordwise surface vortex distribution is applied. Hess 
reports that his latest 3-D program has a parabolic vorticity 
option (analogous to that of Reference 22) and that for the 
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supercritical geometry the parabolic option would have 
substantially reduced the numerical instability. Nonetheless, 
the present example does support the conclusion that combined 
source-doublet methods are significantly more reliable than 
source methods for thin geometries subjected to strong pressure 
loading . 




Figure 51. Supercritical Wing 
Effect of Panel Spacing, Rectangular Planform, AR = 100 

The second example also involves modelling a 2 -d airfoil 
as a high aspect ratio wing. The Williams two-element airfoil 
of Figure 16 was selected in order to determine whether the wake 
from the forward element disturbs the rear element pressure 
distribution. In viscous flow, the geometry of the forward 
element wake is significant; however, for purely inviscid flow 
with weak spanwise gradients, the pressure distribution is anti- 
cipated to be nearly independent of wake shape. The reasonable 
but crude wake geometry of Figure 52 was selected, and the wake 
was allowed to trail several dozen chords downstream. The plan- 
form and spanwise panel spacing definition are similar to those 
of the preceding example. 
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36 PANELS PER STREAMWISE SECTION 





Figure 52. Williams Two-Element Wing Geometry 
Rectangular Planform AR = 100 


For the present method, the chordwise spacing of Figure l6 
was applied. For the Hess program, it was necessary to adjust 
the spacing of the rear element trailing edge panels such that 
the upper and lower surface control points nearly aligned. 

Without this adjustment, the calculated lift coefficient based 
on the Hess program equal pressure Kutta condition is approxi- 
mately 25% too low in comparison with the exact value. The calcu- 
lated chordwise pressure distributions near the root presented in 
Figure 53 show that the present method and exact solutions 
compare well. Even with the adjusted trailing edge panels, the 
Hess program significantly underpredicts lift. The probable 
explanation is that the trailing edge Kutta condition is 
sensitive to the locally strong source strengths. 
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Figure 53. Williams 2-Element Wing Pressure Distribution 


In each of the first two examples, the present method 
generates a slightly lower lift coefficient than the exact 2-D 
solution. This is characteristic of the effects of spanwise 
gradients on a finite wing. For a rectangular wing of aspect 
ratio 100, downwash reduces the root section lift coefficient one 
to two percent below the two-dimensional level. 
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Figures 54-56 present the calculated chordwise pressure 
distributions for an unswept wing of rectangular planform, aspect 
ratio six, and NACA 0012 cross section. Using the symmetry plane 
option, the wing semi-span was panelled, with 8 equal size strips 
spanwise and 40 panels per strip streamwise. Both the Hess and 
present methods show good agreement with experimental data at 
6.75° angle of attack (Reference 24) . No viscous corrections 
have been incorporated. It is noteworthy that in the tip region 
the Hess program results are superior to the present method. In 
particular, the upper and lower surface pressure distribution of 
the present method intersect at approximately 70% chord. It is 
believed that this behavior is associated with the quadratic 
surface fit to the doublet distribution. At the tip strip of 
panels the spanwise distance to the adjacent strip control points 
is twice the distance to the tip edge control points. This non- 
equal spacing, coupled with the large doublet spanwise gradients 
at the tip, generates significant errors in the surface fitting 
algorithm. Future implementation of the improved doublet fit 
approach described in the section "Numerical Solution Formulation" 
is expected to eliminate the difficulty. Of course, the less 
efficient approach of increasing the panel spanwise density near 
the tip would also improve the results. 

The flow around the panelled sphere of Figure 57 was Calcu- 
lated by the present method. Only one-half of the sphere was 
panelled, and the symmetry plane option was employed. A hole 
approximately one percent of local panel dimensions was left 
at the north and south poles to prevent the exact coincidence 
of two edge control points of triangular panels. The flow was 
solved for two free stream directions, one parallel to the x-axis 
and one to the z-axis. The latter solution (Figure 58) is of 
greater interest because non-axisymmetric panelling is used to 
model an axisymmetric flow. For both free stream directions 
there is good agreement with the exact analytical solution. It 
is noteworthy that there is a slight scatter in the calculated 
results for panels at the edge of a section, although the magni- 
tude of the scatter is smaller than the symbol size in Figure 58. 
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Apparently the edge control points do not generate surface 
continuity properties as well as the direct imposition of the 
quadratic surface fit between adjacent panels. The pressure 
distribution for the Hess program was also calculated and agrees 
well with the exact solution, but is not shown here. 

Figure 59 illustrates an isolated axisymmetric fuselage with 
an open base. Identical panelling was used for both the present 
and Hess methods, and the calculated pressure distribution was 
compared to the solution generated by Hess' higher order axisym- 
metric surface singularity program (Reference 25) . The higher 
order solution, considered to be virtually exact, is nearly 
matched by both the present and Hess program solutions. 


AR = 6 X = 1 A = 0 NACA0012 



Root 
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Figure 55. Rectangular Wing Pressure Distribution 

Midspan 

AR = 6 X=1 A = 0 NACA0012 
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Figure 57. Sphere Paneling 
400 Panels per Semisphere 
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Figure 58. Pressure Distribution for a Sphere 


[•^1 







CROSS SECTION 



Figure 59. RML51 F07 Body Alone 

The capability to predict internal flow properties was 
assessed by testing the panelled duct of Figure 60. The geometry 
represents the body of revolution generated by wrapping a NACA 
0010 airfoil around a circular cylinder where the cylinder 
length-to-radius ratio is ten. The minimum internal cross- 
sectional area is a factor of four smaller than the corresponding 
inlet and exit areas. The panel dimensions were chosen to be 
unusually large in comparison to the internal diameter in order 
to amplify numerical difficulties, thereby simplifying the 
assessment of the methods. A cylindrical wake extending 
approximately ten duct lengths downstream was panelled to allow 
Kutta condition satisfaction at the exit. 
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Figure 60. Duct Paneling 

Two Views 


The internal pressiire distribution is shown in Figure 61 for 
the present and Hess programs. To verify that the solution from 
the Hess higher order axisymmetric program is sufficiently close 
to exact to be so designated, both the axial panel distribution 
of Figure 60 and a distribution of double density were tested. 

The present method calculations agree quite well with the exact 
solution. It is noteworthy that the minimum pressure coefficient 
is Cp = -13, which compares to the value Cp = -15 corresponding 
to one-dimensional flow theory applied to a constriction ratio of 
four. The Hess program vastly underpredicts the magnitude of the 
internal pressures. This is characteristic of low order source 
methods, which exhibit significant leakage between control points 
for internal flow. 
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Figure 61. Duct Internal Pressure Distribution 
Exit Kutta Condition Applied 

Because it is conceivable that inaccuracies associated with 
exit Kutta condition satisfaction could have affected the duct 
solution of Figure 61, the example was repeated with no wake 
panelling and no Kutta condition (Figure 62) . Although the result- 
ing magnitude of internal flow velocity is significantly reduced, 
the relative prediction accuracy of the different methods is not 
substantially changed. 

Hess has demonstrated that the inclusion of higher order 
terms significantly improves source method internal flow predic- 
tion accuracy (Reference 26) . Of course, increasing panel density 
also increases the accuracy. For the Hess program, doiobling 
the number of panels per unit length (quadrupling the total number) 
reduces the leakage by an approximate factor of two. It is 
interesting that even without higher order corrections , the 
present method provides satisfactory internal flow solutions for 
low panel density. 
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Figure 62. Duct Internal Pressure Distribution 
Zero Circulation 




The objective of the final example is to determine whether 
the interference effect of a wing on a fuselage pressure distri- 
bution can be predicted without extending the wing panelling 
through the fuselage interior. Whereas the Hess program uses 
internal fuselage panelling, the present method relies entirely 
on surface panels to generate the fuselage lift. A disadvantage 
of the internal panelling approach is that it is not clear 
how various high or low wing configurations should be modelled. 

Figure 63 illustrates the wing-body geometry selected for 
this example. The conventional mid-wing configuration is 
characteristic of geometries for which the Hess program predic- 
tions are very accurate, thereby providing a standard of compari- 
son for the present method. Furthermore, experimental pressure 
data are available (Reference 27) , which provide a second standard. 
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Wing Details 
Airfoil Section 
NACA 65A006 
Aspect Ratio 4 

Taper Ratio 0.6 



The geometry is symmetric with respect to both the wing 
chord plane and the zero butt line plane. Panelling was 
established on only one side of the latter plane of symmetry, 
and an appropriate symmetry option was used in program computa- 
tions. Outboard of the fuselage, the wing was modelled by nine 
streamwise strips of approximately uniform width. Each strip 
contained twenty upper and twenty lower surface panels . The 
fuselage cross section panelling was equally spaced, with each 
panel subtending an arc of 30°. Nineteen fuselage station 
locations were used to define the fuselage panelling. For the 
present method only, the fuselage panel density in the region of 
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the wing intersection was increased and the panel distributions 
modified such that the wing and fuselage panel edges aligned. 
Whether this precaution is beneficial was not determined prior 
to the publication of this report. 

In order to compensate for a deficiency in the present 
program/ it was necessary to model the fuselage afterbody as 
a cylinder. For consistency, the modelling was also used in 
the Hess program. The deficiency is that if the centroids of a 
streamwise strip of wake panels are not colinear, the trailing 
vortex filaments will not be streamwise. In most cases it is 
possible to panel the wake by parallelograms, which eliminates 
the difficulty. Nonetheless, the deficiency in the present 
program should be eliminated. 

The flow around the wing-body was calculated at 4° angle of 
attack. Neither compressibility nor viscous corrections were 
made to the calculated pressures. The test data corresponds to 
free stream Mach number 0.60. 

The calculated and experimental fuselage pressure distribu- 
tions are presented in Figure 64. In spite of significantly 
different approaches used to generate fuselage lift carryover, 
the predictions by the two programs are in close agreement. 

In fact, the total configuration lift coefficients agree to 
three significant figures (0.256 at 4°a). For the Hess pro- 
gram, the slightly milder gradients at the wing root leading 
edge are probably due to the less dense fuselage panelling. The 
fact that the present method pressure distribution is of slightly 
lower magnitude immediately above and below the wing root is 
believed to be attributable to excursions from the nominal 
azimuthal angles 0 of Figure 64 (0 = 45°, 75°, 105°, 135°). 

For the present method, the panelling at the wing-fuselage 
intersection is such that the actual values of 0 are approximate- 
ly 70° and 110° instead of the nominal values 75° and 105°, re- 
spectively. Both the present and Hess methods agree reasonably 
well with the experimental data, in spite of the fact that 
neither viscous nor compressibility corrections were applied. 
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By running the Hess program with and without the fuselage 
afterbody curvature, it was established that more than half of 
the difference between calculated and experimental pressures 
is attributable to modelling the afterbody as a cylinder. 



Fraction Fuselage Length 


Figure 64. Fuselage Pressure Distribution in Presence of Wing 



The wing section pressure distributions at three spanwise 
stations are presented in Figure 65-67. The two methods agree 
well with one another, except near the wing tip. As discussed 
earlier, a more sophisticated doublet surface fitting procedure 
is needed to improve the reliability of the present method at 
wing tips. Over the remainder of the wing, the greater part 
of the differences between calculated and experimental pressures 
is attributable to neglecting the fuselage afterbody curvature. 

The wing-body example demonstrates that the panel edge 
boundary condition control points can provide the proper fuse- 
lage lift carryover without the extension of wing panelling 
through the fuselage. 

Based on the above examples, some deficiencies of the 
present method which can and should be eliminated have been 
uncovered. The only one not discussed above pertains to the 
calculation of the influence coefficients. Only near field 
formulae are currently used. It is believed that this increases 
the time spent computing influence coefficients by an order of 
magnitude over what would be required if both far and intermediate 
formulae were used. As a result, the total computing time for a 
wing-body configuration is approximately five times greater than 
what would be required. This is reflected by comparing the com- 
puting times of the Hess and present programs. For example, on 
the CDC CYBER 175, the duct solution required 100 seconds com- 
puting time in the Hess program versus nearly 400 seconds for 
the present method. In each program, 360 body panels were used 
on each half of the symmetry plane. Clearly, far field formulae 
need to be introduced. 

In general, the above examples confirm that the combined 
source-doublet distribution of Green's identity coupled with 
internal perturbation potential boundary conditions provides 
accurate and stable numerical flow field predictions for a wide 
variety of body shapes, including thin wings and duct interiors. 
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CONCLUSIONS AND RECOMMENDATIONS 

The use of mild combined source-doviblet surface singularity 
distributions eliminates the numerical instabilities associated 
with the excessively strong singularity magnitude of source only 
solutions for highly loaded thin geometries. Also, significantly 
improved prediction accuracy is achieved at sharp concave corners. 
Even without higher order corrections, the use of internal poten- 
tial boundary conditions produces surprisingly good accuracy for 
a wide variety of body geometries. The implementation of Green's 
identity offers several other benefits. Nonlifting configurations 
are made lifting by the simple introduction of a wake. The con- 
figuration doublet distribution automatically readjusts to generate 
lift, thereby eliminating the requirement for separate lifting and 
nonlifting formulations. Internal potential boundary conditions 
can be applied, which tend to improve numerical stability further 
and, as demonstrated in two dimensions, lead to a simplified 
reliable solution method for design problems in which the pressure 
distribution is prescribed. All boundary flow properties can be 
evaluated directly from the local singularity strengths, which 
saves the effort of summing the individual influences of the 
singularities on all panels. This feature also simplifies the 
implementation of mixed Neumann-Dirichlet boundary conditions. 

Two areas of further study are recommended, specifically, 
the improvement of doiiblet continuity properties and the imple- 
mentation of a design formulation for arbitrary bodies. 

It is expected that doublet continuity properties can be 
improved significantly by using common values at adjacent panel 
edges prior to boundary condition satisfaction. At the expense 
of slightly greater computing time, this technique should both 
increase prediction accuracy and reduce the sensitivity to the 
panel distribution. Elimination of the trapezoidal panel restric- 
tion by implementing more comprehensive influence functions should 
also improve continuity characteristics. 

For inverse design problems in which the geometry most 
nearly corresponding to a prescribed pressure distribution is 
to be determined, it is expected that the characteristics of 
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efficiency and stability demonstrated by the two-dimensional 
method will carry over in an application to three-dimensional 
geometries. The two-dimensional method exhibits unusual 
insensitivity to starting geometry and good convergence charac 
teristics even in the difficult leading edge regions. 


McDonnell Aircraft Company 
McDonnell Douglas Corporation 
St. Louis, Missouri 63166 
January 13, 1978 
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